arXiv:2201.07727v1 [astro-ph.EP] 19 Jan 2022 


MNRAS 000, 1-22 (2022) Preprint 20 January 2022 Compiled using MNRAS IATEX style file v3.0 


Investigating the architecture and internal structure of the TOI-561 system 
planets with CHEOPS, HARPS-N and TESS 


G. Lacedelli!2*, T. G. Wilson?, L. Malavolta!2, M. J. Hooton?, A. Collier Cameron?, Y. Alibert^?, 

A. Mortier®’, A. Bonfanti®, R. D. Haywood?, S. Hoyer!®, G. Piotto! 2, A. Bekkelien!!, 

A. M. Vanderburg!*!3, W. Benz^?, X. Dumusque!!, A. Deline!!, M. López-Morales!^, L. Borsato?, 

K. Rice'?/6, L. Fossati?, D. W. Latham!*, A. Brandeker!’, E. Poretti!*:?, S. G. Sousa”, 

A. Sozzetti?! , S. Salmon!!, C. J. Burke"?, V. Van Grootel?, M. M. Fausnaugh!”, V. Adibekyan/?, 

C. X. Huang!??3, H. P. Osborn>-!”, A. J. Mustill^, E. Pallé??, V. Bourrier!!, V. Nascimbeni?, 

R. Alonso???6, G. Anglada??5, T. Bárczy”, D. Barrado y Navascues??, S. C. C. Barros??!, 

W. Baumjohann?, M. Beck!!, T. Beck*, N. Billot!!, X. Bonfils??, C. Broeg??, L. A. Buchhave??, 

J. Cabrera?^, S. Charnoz?, R. Cosentino!®, Sz. Csizmadia?^, M. B. Davies, M. Deleuil!?, L. Delrez2??7, 
O. Demangeon?”?!, B.-O. Demory?, D. Ehrenreich!! , A. Erikson?^, E. Esparza-Borges~>®, 

H.-G. Florén!7?5, A, Fortier*>, M. Fridlund??*°, D. Futyan!!, D. Gandolfi*!, A. Ghedina'?, M. Gillon?, 
M. Giidel*?, P. Gutermann!9^, A, Harutyunyan!?, K. Heng?^^, K. G. Isaak*>, J. M. Jenkins^6, 

L. Kiss^^^5, J. Laskar*?, A. Lecavelier des Etangs??, M. Lendl!!, C. Lovis!!, D. Magrin?, 

L. Marafatto?, A. F. Martinez Fiorenzano!®, P. F. L. Maxted?!, M. Mayor! !, G. Micela??, 

E. Molinari’, F. Murgas”, N. Narita??9?^55, G. Olofsson”, R. Ottensamer?, I. Pagano”®, 

A. Pasetti?, M. Pedani!®, F. A. Pepe!!, G. Peter’, D. F. Phillips!^, D. Pollacco“, 

D. Queloz®!!, R. Ragazzoni!, N. Rando’, F. Ratti??, H. Rauer?^606!. J. Ribas?”28, 

N. C. Santos?9?!, D. Sasselov!^, G. Scandariato>®, S. Seager!?92.99, D. Ségransan!!, L. M. Serrano*!, 
A. E. Simon*, A. M. S. Smith?^, M. Steinberger?, M. Steller’, Gy. Szabó9^65. N. Thomas‘, 

J. D. Twicken 966, S. Udry!!, N. Walton’, J. N. Winn? 


Affiliations are listed at the end of the paper 


Accepted 2022 January 18. Received 2022 January 14; in original form 2021 November 25 


© 2022 The Authors 


2  G.Lacedelli et al. 


ABSTRACT 

We present a precise characterization of the TOI-561 planetary system obtained by combining previously published data with 
TESS and CHEOPS photometry, and a new set of 62 HARPS-N radial velocities (RVs). Our joint analysis confirms the presence 
of four transiting planets, namely TOI-561 b (P = 0.45 d, R = 1.42 Re, M = 2.0 Mẹ), c (P = 10.78 d, R = 2.91 Re, 
M = 5.4 Mo), d (P = 25.7 d, R = 2.82 Rẹ, M = 13.2 Mo) and e (P = 77 d, R = 2.55 Re, M = 12.6 Rẹ). Moreover, we 
identify an additional, long-period signal (> 450 d) in the RVs, which could be due to either an external planetary companion 
or to stellar magnetic activity. The precise masses and radii obtained for the four planets allowed us to conduct interior structure 
and atmospheric escape modelling. TOI-561 b is confirmed to be the lowest density (py = 3.8 + 0.5 g cm?) ultra-short period 
(USP) planet known to date, and the low metallicity of the host star makes it consistent with the general bulk density-stellar 
metallicity trend. According to our interior structure modelling, planet b has basically no gas envelope, and it could host a certain 
amount of water. In contrast, TOI-561 c, d, and e likely retained an H/He envelope, in addition to a possibly large water layer. 
The inferred planetary compositions suggest different atmospheric evolutionary paths, with planets b and c having experienced 
significant gas loss, and planets d and e showing an atmospheric content consistent with the original one. The uniqueness of 
the USP planet, the presence of the long-period planet TOI-561 e, and the complex architecture make this system an appealing 
target for follow-up studies. 


Key words: stars: individual: TOI-561 (TIC 377064495, Gaia EDR3 3850421005290172416) — techniques: photometric — 


techniques: radial velocities — planets and satellites: fundamental parameters — planets and satellites: interiors 


1 INTRODUCTION 


Since the announcement of the first exoplanet orbiting a Sun-like 
star (Mayor & Queloz 1995), the growing number of discoveries in 
exoplanetary science have yielded a surprising variety of exoplan- 
ets and exoplanetary systems. The field has benefited hugely from 
dedicated space-based missions, such as CoRoT, Kepler, K2 (Baglin 
et al. 2006; Borucki et al. 2010; Howell et al. 2014), and recently 
TESS (Ricker et al. 2014). With more than 170 confirmed planets, 
and ~ 4000 planet candidates, the majority of which will likely turn 
out to be planets, TESS has increased the census of confirmed exo- 
planets to more than 4500! . Alongside the aforementioned missions, 
which are designed to discover a large number of exoplanets by 
searching for transit-like signatures around hundreds of thousands of 
stars, new characterization missions, with a specific focus on the de- 
tailed study of known exoplanets, are now starting to operate. Among 
them, the CHaracterising ExOPlanet Satellite (CHEOPS, Benz et al. 
2021), launched on 18 December 2019, is a 30 cm telescope which 
is collecting ultra-high precision photometry of known exoplanets, 
aiming at their precise characterization. CHEOPS met its precision 
requirements both on bright and faint stars, achieving a noise level 
of ~ 15 ppm per 6 h intervals for V ~ 9 mag stars, and 75 ppm per 
3 h for V ~ 12 mag stars (Benz et al. 2021). The importance of such 
a high photometric precision is reflected in CHEOPS' first scientific 
results, which span a variety of different fields (Lendl et al. 2020; 
Bonfanti et al. 2021b; Leleu et al. 2021; Delrez et al. 2021; Morris 
et al. 2021a; Borsato et al. 2021; Van Grootel et al. 2021; Szabó 
et al. 2021; Hooton et al. 2021; Swayne et al. 2021; Maxted et al. 
2021; Barros et al. 2022; Wilson et al. 2022; Deline et al. 2022). 
As part of its main scientific goals, CHEOPS is refining the radii of 
known exoplanets to achieve the precision on the bulk density needed 
for internal structure and atmospheric evolution modelling. To fulfil 
this aim, radial velocity (RV) follow-ups using high-precision spec- 
trographs are essential to provide the precise planetary masses that 
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can be combined with radii measurements to determine accurate 
densities. Among the exoplanets having both radius and mass mea- 
surements, the ones in well-characterised multiplanetary systems are 
of particular interest, since they allow for investigation of their forma- 
tion and evolution processes through comparative planetology, e.g. 
by comparing their individual inner bulk compositions (e.g. Guen- 
ther et al. 2017; Prieto-Arranz et al. 2018), by studying their mutual 
inclinations and eccentricities (e.g. Fabrycky et al. 2014; Van Eylen 
et al. 2019; Mills et al. 2019), and by investigating the correlations 
between their relative sizes, masses and orbital separations (e.g. Lis- 
sauer et al. 2011; Ciardi et al. 2013; Millholland et al. 2017; Weiss 
et al. 2018; Jiang et al. 2020; Adams et al. 2020). 


Within this context, TOI-561, announced simultaneously by 
Lacedelli et al. (2021) and Weiss et al. (2021) (L21 and W21 here- 
after, respectively), is a particularly interesting system, both from 
the stellar (Section 2.1) and planetary (Section 2.2) perspective. The 
low stellar metallicity, the presence of an ultra-short period (USP) 
planet, where USP planets are meant here as planets with periods 
shorter than one day and radii smaller than 2 Rg, and the complex- 
ity of its planetary configuration make TOI-561 an appealing target 
for in-depth investigations. In this study, we combine literature data 
with new TESS observations (Section 3.1), CHEOPS photometry 
(Section 3.2), and HARPS-N RVs (Section 3.3) to shed light on the 
planetary architecture and infer the internal structure of the transiting 
planets. After assessing the planetary configuration using CHEOPS 
observations (Section 4.1) and performing a thorough analysis of 
the global RV data set (Section 4.2), we jointly modelled the pho- 
tometric and spectroscopic data to obtain the planetary parameters 
(Section 5). We used our derived stellar and planetary properties to 
model the internal structures of the transiting planets (Section 6) and 
their atmospheric evolution (Section 7), before discussing our results 
and presenting our conclusions (Section 8). 


The TOI-561 system as seen by CHEOPS, HARPS-N and TESS 


2 THE TOI-561 SYSTEM 


2.1 The host star 


TOI-561 is an old, metal-poor, thick disk star (L21, W21), slightly 
smaller and cooler than the Sun, located ~ 84 pc away from the Solar 
System. We report the main astrophysical properties of the star in 
Table 1. 

We adopted the spectroscopic parameters and stellar abundances 
from L21 (Table 1), which were derived exploiting the high SNR, 
high-resolution HARPS-N co-added spectrum (L21, § 3.1) through 
an accurate analysis using three independent methods, namely the 
ARES+MOOG equivalent width method (Sousa 2014; Mortier et al. 
2014), the Stellar Parameter Classification (Buchhave et al. 2012, 
2014) and the CCFpams method (Malavolta et al. 2017). 

Taking advantage of the updated parameters coming from the 
Gaia EDR3 release (Gaia Collaboration et al. 2021), we then used 
the L21 spectral parameters as priors on spectral energy distribution 
selection to infer the stellar radius (Ry) of TOI-561 using the infrared 
flux method (IRFM, Blackwell & Shallis 1977). The IRFM compares 
optical and infrared broadband fluxes and synthetic photometry of 
stellar atmospheric models, and uses known relationships between 
stellar angular diameter, Te and parallax to derive Rx, in a MCMC 
fashion as detailed in Schanche et al. (2020). For this study, we 
retrieved from the most recent data releases the Gaia G, Ggp, Grp 
(Gaia Collaboration et al. 2021), 2MASS J, H, K (Skrutskie et al. 
2006), and WISE W1, W2 (Wright et al. 2010) broadband photometric 
magnitudes, and we used the stellar atmospheric models from the 
ATLAS Catalogues (Castelli & Kurucz 2003) and the Gaia EDR3 
parallax with the offset of Lindegren et al. (2021) applied, to obtain 
R, = 0.843 + 0.005 Ro. 

We combined two different sets of stellar evolutionary tracks and 
isochrones, PARSEC2 (PAdova & TRieste Stellar Evolutionary Code, 
v1.2S; Marigo et al. 2017) and CLES (Code Liégeois d'Évolution 
Stellaire, Scuflaire et al. 2008), to derive the stellar mass (Mx) and 
age (tx) of TOI-561. As the star is significantly alpha-enhanced, 
we avoided using [Fe/H] as a proxy for the stellar metallicity; in- 
stead, we inserted both [Fe/H] and [o/Fe] in relation (3) provided 
by Yi et al. (2001), obtaining an overall scaling of metal abundances 
[M/H] = -0.23 + 0.06. Besides [M/H], the main input parameters 
for computing M, and t were Teg and Rw. In addition, we used as in- 
puts log Rox and the upper limit on v sin i from L21, and the yttrium 
over magnesium abundance [Y/Mg] = —0.22 + 0.07, as computed 
from [Mg/H] and [Y/H] reported by W21. These indices improve 
the model convergence by discarding unlikely young isochrones, as 
broadly discussed in Section 2.2.3 of Bonfanti & Gillon (2020), 
and references therein. The PARSEC results were obtained using 
the isochrone placement algorithm of Bonfanti et al. (2015, 2016), 
which retrieves the best-fit parameters by interpolating within a pre- 
computed grid of models, while the CLES algorithm models di- 
rectly the star through a Levenberg-Marquardt minimisation (Salmon 
et al. 2021). The final adopted values (M, = 0.806 + 0.036 Mo, 
tk = 110725 Gyr) are a combination of the outputs from both sets 
of models, as described in detail in Bonfanti et al. (2021b). The 
derived mass and radius, listed in Table 1, are consistent within 
lo with the values reported in L21 (Ry = 0.849 + 0.007 Ro, 
My = 0.785 + 0.018 Mo). 


? http://stev.oapd.inaf.it/cgi-bin/emd 


Table 1. Stellar properties of TOI-561. 
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TOI-561 

TIC 377064495 
Gaia EDR3 3850421005290172416 
2MASS J09524454+0612589 
Parameter Value Source 
RA (J2016; hh:mm:ss.ss) 09:52:44.43 A 
Dec (J2016; dd:mm:ss.ss) —-06:12:57.94 A 
Ha (mas yr!) —108.504+0.022 A 
us (mas yr!) -61.279+0.019 A 
y (kms7!) 79.54 + 0.56 A 
Parallax (mas) 11.8342 + 0.0208 A 
Distance (pc) 84.25 + 0.12 B 
TESS (mag) 9.527 + 0.006 C 
G (mag) 10.0181 + 0.0028 A 
Gpp (mag) 10.3945 + 0.0028 A 
Grp (mag) 9.4692 + 0.0038 A 
V (mag) 10.252 + 0.006 C 
B (mag) 10.965 + 0.082 C 
J (mag) 8.879 x 0.020 D 
H (mag) 8.504 + 0.055 D 
K (mag) 8.394 x 0.019 D 
WI (mag) 8.337 x 0.023 E 
W2 (mag) 8.396 + 0.020 E 
Teg (K) 5372 + 70 F 
log g (cgs) 4.50 + 0.12 F 

[Fe/H] (dex) —0.40 + 0.05 F 

[Mg/H] (dex) —0.17 x 0.05 F 

[Si/H] (dex) —0.22 x 0.05 F 

[Ti/H] (dex) —0.12 x 0.03 F 

[a/Fe] (dex) 0.23 + 0.04 F 

[M/H] (dex) —0.23 + 0.06 G 

Y/Mg] (dex) —0.22 x 0.07 G^ 
log Rix —5.003 + 0.012 F 
v sini (kms!) 29 F 
Rx (Ro) 0.843 + 0.005 G, IRFM 
My (Mo) 0.806 + 0.036 G, isochrones 
tx (Gyr) 11.0725 G, isochrones 
Px (Po) 1.34 + 0.06 G, from R, and M, 
px (g cm?) 1.89 + 0.09 G, from Rẹ and M, 
Ly (Lo) 0.533 + 0.029 G, from Ry and Ter 
Spectral type G9V F 


A) Gaia EDR3 (Gaia Collaboration et al. 2021). B) Bailer-Jones et al. 
(2021). C) TESS Input Catalogue Version 8 (TICv8, Stassun et al. 2018). 
D) Two Micron All Sky Survey (2MASS, Cutri et al. 2003). E) Wide-field 
Infrared Survey Explorer (WISE; Wright et al. 2010). F) L21. G) This 


work. 


@ Based on W21 abundances. 
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2.2 The planetary system 


The discovery of a multiplanetary system orbiting TOI-561 was an- 
nounced simultaneously by L21 and W21 in two independent papers. 
The main planetary parameters from both studies are reported in Ta- 
ble 2. 

The two papers presented different RV data sets, collected with 
HARPS-N and HIRES respectively, to confirm the planetary nature 
of three candidates identified by TESS in sector 8, the only available 
sector at the time of the publications. The TESS-identified signals had 
periods of ~ 0.45, ~ 10.8, and ~ 16 days. The two inner candidates 
were confirmed by both L21 and W21, with the names of TOI-561 
b (an USP super-Earth, with period Py ~ 0.4465 d, and radius Rp ~ 
1.4 Ra), and TOI-561 c (a warm mini-Neptune, with Pe ~ 10.779 d, 
and Re ~ 2.9 Ro). However, two different interpretations for the 
third TESS signal were proposed by the authors. In the scenario 
presented in L21, the two transits related to the third TESS signal 
were interpreted as single transits of two distinct planets, TOI-561 
d (Pq ~ 25.6 d, Ra ~ 2.5 Re), and TOI-561 e (Pe ~ 77 d, Re ~ 
2.7 Re). The periods of these two planets were inferred from the 
RV analysis, which played an essential role in determining the final 
planetary architecture. In fact, the ephemeris match between the 
RV and photometric fits (See Fig. 5 of L21) and the non-detection 
of the 16 d signal in the HARPS-N data set, combined with the 
different durations of the two TESS transits and results from the 
long-term stability analysis led the authors to converge on a 4-planet 
configuration, presenting robust mass and radius detection for all the 
four planets in the system (L21, Table 5). In contrast, W21 proposed 
the presence of a single planet at the period suggested by TESS (TOI- 
561 f, Pp ~ 16.29 d, Rf ~ 2.3 Ro), based on the analysis of the two 
available transits. W21 pointed out that the 8.1 d alias of planet f’s 
orbital period is also consistent with the TESS data, with the even 
transit falling into the TESS download gap, even though in this case 
the transit duration would be too long compared to what is expected 
for a 8 d period planet ($4.9, W21). However, the authors could 
not obtain an accurate mass determination for this planet, with the 
60 HIRES RVs being consistent with a non-detection (8 7.2, W21). 
An additional discrepancy between the two studies is the mass of 
the USP planet, differing by almost a factor two. According to the 
W21 analysis, TOI-561 b has a mass of 3.2 + 0.8 Mẹ, making 
it consistent with a rocky composition and placing it among the 
population of typical small (< 2 Rẹ), extremely irradiated USP 
planets (Sanchis-Ojeda et al. 2015; Dai et al. 2021). Instead, assuming 
the low mass (My = 1.59 + 0.36 Ma) inferred from L21 analysis, 
TOI-561 b is not consistent with a pure rocky composition, and it 
is the lowest density USP super-Earth known to date, calling for 
a more complex interpretation (e.g. lighter core composition, deep 
water reservoirs, presence of a high-metallicity, volatile materials or 
water steam envelope, etc.). 

The complexity of this system and the differences between the two 
studies demanded further investigations. We therefore decided to col- 
lect additional, precise photometric and RV data (Section 3) to shed 
light on the planetary configuration and on the internal composition 
of the TOI-561 planets. 


3 OBSERVATIONS 
3.1 TESS photometry 


During its two-year primary mission (Ricker et al. 2014), TESS ob- 
served TOI-561 in two-minute cadence mode between 2 February 
and 27 February 2019 (sector 8). After entering its extended mission, 
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Table 2. Literature parameters of the proposed planets orbiting TOI-561. 


TOI-561 b Lacedelli et al. (2021) Weiss et al. (2021) 
P (d) 0.446578 + 0.000017 — 0.446573*0 00009 
To (TBJD) 1517.498 + 0.001 1517.4973 + 0.0018 
Rp (Re) 1.423 + 0.066 1.45 x 0.11 

K (ms?!) 1.56 + 0.35 3.12 0.8 

My (Me) 1.59 + 0.36 3.2 € 0.8 

TOI-561 c 

P (d) 10.779 + 0.004 10.77892 + 0.00015 
To (TBJD) 1527.060 + 0.004 1527.05825 0.00053 
Rp (Re) 2.878 + 0.096 2.90 + 0.13 

K (ms?!) 1.84 + 0.33 2.4 x 0.8 

My (Me) 5.40 + 0.98 7.0 € 2.3 

TOI-561 d 

P (d) 25.62 + 0.04 : 

To (TBJD) 1521.882 + 0.004 : 

Rp (Re) 2.53 + 0.13 : 

K (ms?!) 3.06 + 0.33 - 

Mp (Mo) 11.95 + 1.28 - 

TOI-561 e 

P (d) 71.23 + 0.39 - 

To (TBJD) 8538.181 + 0.004 - 

Rp (Re) 2.67+0.11 - 

K (ms!) 2.84 + 0.41 - 

Mp (Mo) 16.0 + 2.3 : 

TOI-561 f 

P (d) - 16.287 + 0.005 

To (TBID) | - 1521.8828 + 0.0035 
Rp (Re) - 2.32 + 0.16 

K (ms?!) - 0.9 + 0.6 

Mp Me) — - 3.085 

M. (Mo) 0.785 + 0.018 0.805 + 0.030 


a Referred as TOI-561 d in W21. 


TESS re-observed the star in two-minute cadence mode during sector 
35, between 9 February and 6 March 2021. At the beginning of the 
second orbit, the spacecraft dropped out of Fine Pointing mode for 
3.44 days, entering Coarse Pointing mode?. Data collected during 
Coarse Pointing mode were flagged and removed from the Pre-search 
Data Conditioning Simple Aperture Photometry (PDCSAP, Smith 
et al. 2012; Stumpe et al. 2012, 2014) light curves, leading to a total of 
19.86 days of science data. The photometric observations of TOI-561 
were reduced by the Science Processing Operations Center (SPOC) 
pipeline and searched for evidence of transiting planets (Jenkins et al. 
2016; Jenkins 2020). For our photometric analysis, we used the light 
curves based on the PDCSAP, downloading the two-minute cadence 


3 See TESS Data Release Notes: Sector 35, DR51 (https://archive. 
stsci.edu/tess/tess_drn.html). 
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Table 3. Number of TOI-561 transits observed by TESS. 


TOI-561b  TOILS6Ic  TOI-561d TOI-561 e 
Sector 8 41 2 1 1 
Sector 35 43 1 - - 


data from the Mikulski Archive for Space Telescopes (MAST)*, and 
removing all the observations encoded as NaN or flagged as bad- 
quality (DQUALITY>9) points by the SPOC pipeline’. We performed 
outlier rejection by doing a cut at 30 for positive outliers and 5c 
(i.e. larger than the deepest transit) for negative outliers. The result- 
ing TESS light curves of sectors 8 and 35 are shown in Figure 1, and 
Table 3 summarizes the total number of transits observed by TESS 
for each planet. 

To refine the ephemeris of planet d in time for the scheduling 
of the CHEOPS observations (Section 4.1), we also extracted the 
10-minute cadence light curve of sector 35 using the quick-look 
TESS Full Frame Images (FFIs) calibrated using the TESS Image 
CAlibrator® package (tica, Fausnaugh et al. 2020). 


32 CHEOPS photometry 


To confirm the planetary architecture and improve the planetary pa- 
rameters, we obtained three visits of TOI-561 with CHEOPS, the 
ESA small class mission dedicated to the characterization of known 
exoplanets (Benz et al. 2021). The observations, collected within the 
Guaranteed Time Observing (GTO) programme, were carried out 
between 23 January and 15 April 2021, for a total of 73.85 hours on 
target. During the three visits, we observed a total of eight transits 
of TOI-561 b, two transits of TOI-561 c, and one transit of TOI-561 
d. The three CHEOPS light curves have an observing efficiency, i.e. 
the actual time spent observing the target with respect to the total 
visit duration, of 64%, 75%, and 61%, respectively. The observing 
efficiency is linked to data gaps, which are intrinsically present in 
all CHEOPS light curves (see e.g. Delrez et al. 2021, Bonfanti et al. 
2021b, Leleu et al. 2021), and are related to CHEOPS’s low-Earth 
orbit. In fact, during (1) South Atlantic Anomaly (SAA) crossing, 
(2) target occultation by the Earth, and (3) too high stray light con- 
tamination, no data are downlinked. This results in data gaps, whose 
number and extension depend on the target sky position (Benz et al. 
2021). For all the visits, we adopted an exposure time of 60 s. The 
summary log of the CHEOPS observations is reported in Table 4. 
Data were reduced using the latest version of the CHEOPS auto- 
matic Data Reduction Pipeline (DRP v13; Hoyer et al. 2020), which 
performs aperture photometry of the target after calibrating the raw 
images (event flagging, bias, gain, non-linearity, dark current, and 
flat field) and correcting them for instrumental and environmental 
effects (smearing trails, cosmic rays, de-pointing, stray light, and 
background). The target flux is obtained for a set of three fixed- 
radius apertures, namely R = 22.5 arcsec (RINF), 25.0 arcsec (DE- 
FAULT), 30.0 arcsec (RSUP), plus an additional one specifically 
computed to optimize the radius based on the instrumental noise 
and contamination level of each target (OPTIMAL). Moreover, the 
DRP estimates the contamination in the photometric aperture due to 


4 https://mast.stsci.edu/portal/Mashup/Clients/Mast/ 
Portal .html 

5 https://archive.stsci.edu/missions/tess/doc/ 

EXP- TESS- ARC- ICD- TM-0014-Rev-F.pdf 

6 https://archive.stsci.edu/hlsp/tica 


nearby targets using the sources listed in the Gaia DR2 catalog (Gaia 
Collaboration et al. 2018) to simulate the CHEOPS Field-of-View 
(FoV) of the target, as described in detail in Hoyer et al. (2020). No 
strong contaminants are present in the TOI-561 FoV, and the main 
contribution to the contamination is due to the smearing trails of a 
G = 10.20 mag star at a projected sky distance of ~ 117.9 arcsec, 
which rotates around the target inside the CCD window because of 
the CHEOPS field rotation (Benz et al. 2021). During the third visit 
three telegraphic pixels (pixels with a non-stable and abnormal be- 
haviour during the visit) appeared within the CHEOPS aperture, one 
of them inside the CHEOPS PSF (Figure 2). A careful treatment, 
described in detail in Appendix A, was applied to correct for their 
effect. In the subsequent analysis we adopted for all the visits the 
RINF photometry (see Figure A1 in Appendix A), which minimized 
the light curve root mean square (RMS) dispersion, and we removed 
the outliers by applying a 4c clipping. 

Finally, a variety of non-astrophysical sources, such as varying 
background, nearby contaminants or others, can produce short-term 
photometric trends in the CHEOPS light curves on the timescale 
of one orbit, due to the rotation of the CHEOPS FoV around the 
target and due to the nature of the spacecraft orbit. To correct for 
these effects, we detrended the light curves using the basis vectors 
provided by the DRP, as detailed in Section 5. The resulting detrended 
light curves are shown in Figure 3. 


3.3 HARPS-N spectroscopy 


In addition to the 82 RVs published in L21, we collected 62 high- 
resolution spectra using HARPS-N at the Telescopio Nazionale 
Galileo (TNG), in La Palma (Cosentino et al. 2012, 2014). These were 
used to refine the planetary masses and confirm the system configu- 
ration. The new observations were collected between 15 November 
2020 and 1 June 2021. Following the same strategy of the previous 
season (L21), in addition to 30 single observations, we collected six 
points per night on 8 and 10 February 2021, and two points per night 
on ten additional nights, specifically targeting the USP planet. The 
exposure time for all the observations was set to 1800 s, resulting 
in a signal-to-noise ratio (SNR) at 550 nm of 83 + 20 (median + 
standard deviation) and a radial velocity measurement uncertainty of 
1.0: 0.4 m s7}. All the observations were gathered with the second 
HARPS-N fibre illuminated by the Fabry-Perot calibration lamp to 
correct for the instrumental RV drift. 

We reduced the global HARPS-N data set (144 RVs in total) using 
the new version of the HARPS-N Data Reduction Software based 
on the ESPRESSO pipeline (DRS, version 2.3.1; see Dumusque 
et al. 2021 for more details). We used a G2 flux template to correct 
for variations in the flux distribution as a function of wavelength, 
and a G2 binary mask to compute the cross-correlation function 
(CCF, Baranne et al. 1996; Pepe et al. 2002). We report the RVs 
and the associated activity indices (see Section 4.2) with their 10 
uncertainties in Table 5. As in L21, we removed from the first season 
data set five RVs with associated errors > 2.5 m s^! from spectra 
with SNR « 35 (see Appendix B in L21). All the RV uncertainties of 
the second season data set were below 2.5 m s^, so no points were 
removed. 


3.4 HIRES spectroscopy 


We included in our analysis 60 high-resolution spectra collected 
with the W.M. Keck Observatory HIRES instrument on Mauna Kea, 
Hawaii between May 2019 and October 2020. The data set was 
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Figure 1. TESS sector 8 (left) and 35 (right) PDCSAP light curves of TOI-561. In the top panel, the dark red solid line shows the best-fitting transit and Matérn-3/2 
kernel Gaussian Process (GP) model, as detailed in Section 5. The central panel shows the flattened light curve after the removal of the GP component, with the 
best-fitting transit model superimposed (dark red solid line). The transits of planets c, d and e are labelled and highlighted with orange, green and red vertical 
lines, respectively. The expected locations of the transits of planets c and d occurring during the data gaps of sector 35 are marked with pale, dashed orange and 


green lines, respectively. Planet e is not expected to transit in sector 35. The transits of the USP planet are too shallow to be individually visible, and are not 
indicated. Light curve residuals are shown in the bottom panel. 


Table 4. Log of TOI-561 CHEOPS observations. 


Visit File key Starting date Duration Datapoints Efficiency Exposure time Planets 
(#) (UTC) (h) (#) (%) (s) 
1 CH PRI100031 TG037001 V0200  2021-01-23T15:29:07 15.67 604 64 60 b,c 
2 CH PRI00008 TGO000811 V0200  2021-03-29T10:19:08 4.42 207 75 60 b,c 
3 CH PRI100031 TG039301 V0200  2021-04-12T23:52:28 53.76 1978 61 60 b,d 
Table 5. HARPS-N radial velocity and activity indices measurements. 
BJDTDB RV ORV FWHM OFWHM BIS OBIS Contrast  Ocontr — S-index os Ha CHa 
(d) (m s7!) (ms) (kms!) (kms!) (ms!) (ms) (dex) (dex) 
2458804.70780 79695.97 1.13 6.415 0.002 -86.82 2.26 59.879 0.021 0.1643 0.0005 0.2101 0.0002 
2458805.77552 79699.66 0.85 6.419 0.002 -85.13 1.71 59.810 0.016 0.1702 0.0003 0.2124 0.0001 
2458806.76769 79697.50 0.91 6.415 0.002 -83.66 1.82 59.861 0.017 0.1689 0.0003 0.2082 0.0001 


This table is available in its entirety in machine-readable form. 
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Figure 2. Extraction of 60 x 60 arcsec of the CHEOPS FoV during the third 
visit centered on TOI-561. The dashed black circle represents the RINF pho- 
tometric aperture surrounding the CHEOPS PSF, whose centroid is marked 
by the black cross. The positions of the three identified telegraphic pixels, 
including the one located within the CHEOPS PSF (see Appendix A), are 
highlighted by the red, circled crosses. 


published in W21, and we refer to that paper for details regarding 
the observing and data reduction procedures. The HIRES data set 
has an RMS of 5 m s^! , and a median individual RV uncertainty of 
1.4m s^! (W21). 


4 PROBING THE SYSTEM ARCHITECTURE 
4.1 CHEOPS confirmation of TOI-561 d 


To solve the discrepancy among the planetary architectures proposed 
by L21 and W21 (Section 2.2), we initially looked for the transits 
of TOI-561 d (~ 25 d) and TOI-561 f (~ 16 d) in the TESS sector 
35 light curve, whereas TOI-561 e (~ 77 d) was not expected to 
transit during those TESS observations. However, as shown in the 
top panel of Figure 4, the transits of planet d and f occurred during 
the light curve gap (Section 3.1), and so we could not use the new 
TESS data to conclusively discriminate between the two planetary 
configurations. Nonetheless, these observations ruled out the planet 
f alias at 8.1 d mentioned in W21, since no transit events were 
detected at its predicted transit times. 

We therefore decided to probe the L21 scenario collecting a transit 
of TOI-561 d using CHEOPS. We opted for the scheduling of the 
last seasonal observing window, in April, in order to take advantage 
of the most updated ephemeris to optimize the scheduling. For this 
reason, we performed a global fit adding to the literature data a 
partial set of the new HARPS-N RVs, as of 16 March 2021, and 
including the TESS sector 35 light curve extracted from the second 
data release of the tica FFIs in March 2021. Even if no transit was 
detected, the new TESS sector helped to reduce the time window 
in which to search. In fact, the TESS data partially covered the 30- 
uncertainty transit window, enabling us to exclude some time-spans 
in the computation of the CHEOPS visit. Thanks to the ephemeris 
update, the CHEOPS 30 observing window shrank from ~ 7.4 d 
to ~ 2.2 d, demonstrating the importance of the early TESS data 
releases in the scheduling of follow-up observations. The bottom 


panel of Figure 4 shows the CHEOPS visit scheduled to observe 
TOI-561 d, whose transit occurred almost exactly at the predicted 
time, so confirming the planetary period and giving further credence 
to the 4-planet scenario proposed by L21 . 

Even updating the ephemeris using the partial new HARPS-N data 
set, the last possible CHEOPS observing window of TOI-561 e in the 
2021 season was still longer than seven days because of ephemeris 
uncertainties. Even including the full set of RVs would have not 
helped as the target was no longer observable with CHEOPS when 
the HARPS-N campaign finished. Given the high pressure on the 
CHEOPS schedule, we therefore plan the TOI-561 e observations 
for the 2022 observing season. The ephemeris for the 2022 CHEOPS 
observations will be updated using the TESS Sectors 45 and 46 
observations in Nov-Dec 2021, and the results will be presented in a 
future publication. 


4.2 Additional signals in the RV data 


Before proceeding with the global modelling, we analyzed the RV 
data sets in order to confirm the robustness of the L21 scenario 
and search for potential new signals. The £ |-periodogram" (Hara 
et al. 2017) of the combined HARPS-N and HIRES RVs (Figure 5) 
shows four significant peaks corresponding to the planetary periods 
reported in L21, plus hints of a possible longer period signal with 
a broad peak around 400 — 600 days. We investigated the presence 
of this additional signal in a Bayesian framework using PyORBIT? 
(Malavolta et al. 2016, 2018), a package for light curve and RVs anal- 
ysis. We employed the dynesty nested-sampling algorithm (Skilling 
2004; Skilling 2006; Speagle 2020), assuming 1000 live points, and 
including offset and jitter terms for each data set. We first performed a 
4-planet fit of the combined data sets, using the L21 values to impose 
Gaussian priors on periods and transit times,’ and assuming eccen- 
tric orbits with a half-Gaussian zero-mean prior on the eccentricity 
(with variance 0.098; Van Eylen et al. 2019), except for the circular 
orbit of the USP planet. We let the semi-amplitude K vary between 
0.01 and 100 m s^!. As can be seen in Figure 6, the RV residuals 
show an anomalous positive variation at ~ 9000 BJD-2450000, and 
the Generalized Lomb-Scargle (GLS, Zechmeister & Kürster 2009) 
periodogram of the RV residuals revealed the presence of a signifi- 
cant, broad peak at low frequencies. Moreover, the HARPS-N jitter 
was 1.84 m s~!, which is unusually high when compared to the value 
reported in L21 (cHARPS-N = 1.29 +0.23 m s-1). We therefore per- 
formed a second fit including a fifth Keplerian signal, allowing the 
period to span between 2 and 900 d. According to the Bayesian Ev- 
idence, this model is strongly favoured with respect to the 4-planet 
model, with a difference in the logarithmic evidences Aln Z = 19.0 
(Kass & Raftery 1995).!° Moreover, the HARPS-N jitter decreased 
to ~ 1.37 m s™!. After this fit, the periodogram of the residuals did 
not show evidence of additional significant peaks (Figure 6). This is 
confirmed also by the comparison with a 6-Keplerian model that we 
tested, with the period of the sixth Keplerian free to span between 2 
and 900 d, whose Bayesian Evidence differed by less than 2 from the 
5-Keplerian model one, indicating that there was no strong evidence 
to favour a more complex model (Kass & Raftery 1995). 


7 https://github.com/nathanchara/llperiodogram. 

8 https://github.com/LucaMalavolta/PyORBIT, V8.1. 

? We note that we obtained the same results when using uniform, uninfor- 
mative priors, also for the 5- and 6-Keplerian fits. 

10 According to Kass & Raftery (1995), a difference A In Z > 5 sets a strong 
evidence against the null hypothesis, which in our case corresponds to the 
4-planet model. 
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Figure 3. CHEOPS detrended light curves of TOI-561. Visits 1, 2 and 3 are shown in the top left, top right, and bottom panel, respectively. The best-fitting 
model is over-plotted as a red solid line, and residuals are shown for each visit. The transits of planets b, c, and d are highlighted with blue, orange, and green 


triangles, respectively. 


The fitted period of the fifth Keplerian was ~ 480 d. Such a long- 
term signal could be induced either by stellar activity, considering 
that stellar magnetic fields related to magnetic cycles can show vari- 
ability on timescales of the order of 1 —3 years (e.g. Collier Cameron 
2018; Hatzes 2019; Crass et al. 2021), or by an additional long-period 
planet. We refer here to an eventual long-period planet because, given 
the inferred semi-amplitude of ~ 2m s7! (Table 6), an external com- 
panion with mass equal to 13 Mj (assuming this value to be the 
threshold between planetary and sub-stellar objects) would have an 
inclination of ~ 0.03 deg. Such an inclination would imply an almost 
perpendicular orbit with respect to the orbital plane of the four inner 
planets, hinting at a very unlikely configuration. Therefore, in the 
hypothesis of the presence of an external companion, it would most 
likely be a planetary-mass object. 

On one hand, all the five signals, including the long-term one, are 
recovered in an independent analysis that we performed with the 
CCF-based scALPzrs algorithm (Collier Cameron et al. 2021). Con- 
cisely, SCALPELS projects the RV time series onto the highest variance 
principal components of the time series of autocorrelation functions 
of the CCF, with the aim of distinguishing RV variations caused 
by orbiting planets from activity-induced distortions on each CCF. 
The absence of the signal in the scArPELs shape-driven velocities 
indicates that the long-term periodicity is not due to shape changes 
in the line profiles, supporting the idea of a planetary origin. More- 
over, TOI-561 is not expected to be a particularly active star given 


its old age and low log Rak: as assessed in the L21 and W21 activ- 


MNRAS 000, 1-22 (2022) 


ity analyses. As can be seen in Figure 7, the GLS periodogram of 
the majority of the activity indicators extracted with the HARPS-N 
DRS, i.e. full width at half maximum (FWHM), bisector span (BIS), 
contrast and Hg, do not show significant peaks, with none of them 
exceeding the 0.1 False Alarm Probability (FAP) threshold, which 
we computed using a bootstrap approach, at the frequency of inter- 
est. On the other hand, the periodogram of the S-index, which is 
particularly sensitive to magnetically-induced activity, shows a sig- 
nificant, broad peak at low frequencies, potentially suggesting that 
the previously identified long-term variability is related to stellar 
activity. Considering this, we performed an additional dynesty fit 
assuming a 4-planet model and including a Gaussian Process (GP) 
regression with a quasi-periodic kernel, as formulated in Grunblatt 
et al. (2015), to account for the long-term signal. We modelled si- 
multaneously the RVs and the S-index time series in order to better 
inform the GP (Langellier et al. 2021; Osborn et al. 2021), using two 
independent covariance matrices for each dataset with common GP 
hyper-parameters except for the amplitude of the covariance matrix, 
assuming uniform, non-informative priors on all of them. The fit 
suggests a periodicity longer than ~ 570 d, but the GP model is too 
flexible to derive a precise period value, considering also that the 
global RV baseline (~ 768 d) is comparable with the periodicity of 
the long-term signal. The inferred semi-amplitudes of the four known 
planets differed by less than 0.070 from the 5-Keplerian model ones, 
indicating that the different modelling of the long-term signal is not 
influencing the results for the known, transiting planets. Finally, as 
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Figure 4. Top: 2-min cadence detrended TESS light curve of sector 35. The 
predicted transit times of TOI-561 c and d (according to L21 ephemeris), 
and TOI-561 f (according to W21 ephemeris), are highlighted with orange, 
green and black vertical solid lines, respectively. The black dashed lines 
indicate the predicted position of planet f alias at ~ 8.1 d. The only transit 
present in the light curve is the one of TOI-561 c at ~ 9260 BJD—2450000, 
while no transit events occurred at the predicted times of planet f alias. The 
transits of planet d and f fall into the time series gap. Bottom: CHEOPS visit 
scheduled to observe TOI-561 d. The green vertical solid line indicates the 
predicted transit time used to compute the CHEOPS observing window after 
the ephemeris update (Section 4.1). The transit occurred within the 68 per 
cent highest probability density interval, highlighted by the pale green region. 
We note that this transit is not consistent with the ephemeris propagation of 
planet f, which would have transited at 9319.94 BJD—2450000, so almost 
one day after the observed CHEOPS transit. 


in the case of the 5-Keplerian fit, the HARPS-N jitter is significantly 
improved (cuARPS-N ~ 1.30 m s71) when including the GP model. 
Therefore, since our Bayesian analyses showed that the modelling 
of the long-term signal is necessary to obtain the best picture of the 
system, we decided to perform the global fit assuming a 5-Keplerian 
model, but without drawing conclusions on the origin of the fifth 
signal. We stress that the 5-Keplerian fit does not provide absolute 
evidence of the presence of a fifth planet, since also poorly sampled 
stellar activity could be well modelled using a Keplerian (Pepe et al. 
2013; Mortier & Collier Cameron 2017; Affer et al. 2016), especially 
in our case where the RV baseline is of the order of the signal peri- 
odicity. Since it is not possible to distinguish a true planetary signal 
from an activity signal that has not been observed long enough to 
exhibit a loss of coherence in its phase or amplitude, only a follow-up 
campaign over several years can allow one to better understand the 
nature of this long-term signal. 


5 JOINT PHOTOMETRIC AND RV ANALYSIS 


To infer the properties of the TOI-561 planets, we jointly modelled 
all photometric and spectroscopic data with PyORBIT, using PyDE!! 


1] nttps://github.com/hpparvi/PyDE. 
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Figure 5. £;-periodogram of the combined HARPS-N and HIRES data sets, 
computed on a grid of frequencies from 0 to 2.5 cycles per day. The total 
time-span of the observations is 768 days. The code automatically accounts 
for the offset between HARPS-N and HIRES data by using the mean value 
of each data set. 
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Figure 6. Time series and GLS periodogram of the RV residuals after the 
4-planet and 5-Keplerian fits as described in Section 4.2. In the residuals 
plot, the HARPS-N and HIRES RVs are plotted with dark blue diamonds 
and light blue circles, respectively. In the periodogram plots, the dashed and 
dotted horizontal lines show the 1 and 0.1 per cent False Alarm Probability 
(FAP) level, respectively. The red vertical line indicates the main peak of each 
periodogram. The long-period peak around frequencies 0.0017 — 0.0025 d7! 
(P = 400 — 600 d) in the 4-planet residuals periodogram is modelled by the 
fifth Keplerian, and no more significant peaks are identified in the 5-Keplerian 
residuals periodogram. Moreover, the positive variation at ~ 9000 BJD- 
2450000 in the 4-planet fit residuals disappears in the 5-Keplerian fit residuals. 
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Figure 7. GLS periodogram of the HARPS-N activity indices. The main peak 
of each periodogram is highlighted with a red vertical line. The dashed and 
dotted horizontal lines indicate the 1 and 0.1 per cent FAP levels, respectively. 
The only peak above the 0.1 FAP level is the low-frequency peak in the S- 
index periodogram, as discussed in Section 4.2. 


+ emcee (Foreman-Mackey et al. 2013) as described in Section 5 of 
L21, and adopting the same convergence criteria. We ran 96 chains 
(twice the number of the model parameters) for 250000 steps, dis- 
carding the first 50000 as burn-in. 

Based on the analysis presented in the previous section, we as- 
sumed a 5-Keplerian model, including four planets plus a fifth Ke- 
plerian with period free to span between 2 to 900 d. We fitted a 
common value for the stellar density, using the value reported in 
Table 1 as Gaussian prior. We adopted the quadratic limb-darkening 
law as parametrized by Kipping (2013), putting Gaussian priors on 
the u1, u2 coefficients, obtained for the TESS and CHEOPS passband 
through a bilinear interpolation of limb darkening profiles by Claret 
(2017) and Claret (2021) respectively, and assuming a 1c uncertainty 
of 0.1 for each coefficient. We imposed a half-Gaussian zero-mean 
prior (Van Eylen et al. 2019) on the planet eccentricities, except for 
the USP planet, whose eccentricity was fixed to zero. We assumed 
uniform priors for all the other parameters. 

To model the long-term correlated noise in the TESS light curve, we 
included in the fit a GP regression with a Matérn-3/2 kernel against 
time, as shown in Figure 1, and we added a jitter term to account 
for possible extra white noise. We pre-decorrelated the CHEOPS 
light curves with the pycheops!? package (Maxted et al. 2021), se- 
lecting the detrending parameters according to the Bayes factor to 
obtain the best correlated noise model for each visit. For all the three 
CHEOPS visits, a decorrelation for the first three harmonics of the 


1? https://github.com/pmaxted/pycheops. 
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roll angle was necessary, plus first-order polynomials in time, x-y 
centroid position, and smearing. We then used the detrended light 
curves (Figure 3) for the global PyORBIT fit. In order to check if 
the detrending was affecting our results for the planetary parameters, 
we also performed an independent global analysis with the juliet 
package (Espinoza et al. 2019), including in the global modelling the 
basis vectors selected with pycheops to detrend the data simultane- 
ously. All the results were consistent within lo, indicating that the 
pre-detrending did not significantly alter our inferred results. Finally, 
for both the HARPS-N and HIRES data sets we included jitter and 
offset terms as free parameters. 

We summarize our best-fitting model results in Table 6, and we 
show the transit model, phase-folded RVs and global RV model 
in Figures 8, 9, and 10, respectively. We inferred precise masses 
and radii for all the four planets in the system, whose positions in 
the mass-radius diagram are shown in Figure 11. With a radius of 
Ry = 1.425 + 0.037 Rẹ and a mass of My = 2.00 + 0.23 Me (from 
Ky = 1.93+0.21 m s~!), TOI-561 b is located in a region of the mass- 
radius diagram which is not consistent with a pure rocky composition, 
as will be also shown in Section 6 by our internal structure modelling. 
Our analysis confirms TOI-561 b to be the lowest density (pp = 
3.8+0.5¢ cm?) USP planet known to date (Figure 12). In order to 
further confirm the planetary density, we also performed a specific 
RV analysis of TOI-561 b using the Floating Chunk Offset method 
(FCO; Hatzes 2014). The FCO analysis, detailed in Appendix B, 
confirms the low mass inferred for TOI-561 b, and consequently its 
low density. Thanks to the CHEOPS observations, we also improved 
significantly the radius of TOI-561 c, for which we obtained a value of 


Rc = 2.91+0.04 Re. From the semi-amplitude Ko = 1.81*022 m s7! 


we inferred a mass of Mc = 5390.0 Mə, implying a density of 


pc = 1.2 + 0.2 g cm. From the combined fit of one TESS and one 
CHEOPS transit, we inferred a radius of 2.82 + 0.07 Rg for planet d, 


which has a mass of M; = EE ipe Me (from Kg = gts m s7!) 


and a resulting density of pg = 3.2+0.3 g cm ?. Finally, for TOI-561 
e, which shows a single transit in TESS sector 8, we derived a radius 


of Re = 29578 E Ro, amass of M, = 12.6+ 1.4 Me, and an average 


density of pg = 4.2 + 0.8 g cm. Lastly, the period inferred for the 
fifth Keplerian in the model was 4713+36 d, with a 7.20 detected 


semi-amplitude of 1.94 + 0.27 m sl. As discussed in Section 4.2, 
additional data spanning a longer baseline are needed to definitively 
confirm the planetary nature of this long-term signal. 


6 INTERNAL STRUCTURE MODELLING 


We modelled the internal planetary structure in a Bayesian frame- 
work, following the procedure detailed in Leleu et al. (2021). Our 
model assumes fully-differentiated planets composed of four layers, 
comprising an iron and sulfur central core, a silicate mantle which 
includes Si, Mg and Fe, a water layer, and a pure H/He gas layer. The 
inner core is modelled assuming the Hakim et al. (2018) equation 
of state (EOS), the silicate mantle uses the Sotin et al. (2007) EOS, 
and the water layer uses the Haldemann et al. (2020) EOS. The core, 
mantle and water layer compose the ‘solid’ part of the planet. The 
thickness of the gas envelope is computed as a function of stellar 
age and irradiation, and mass and radius of the solid part, according 
to the model presented in Lopez & Fortney (2014). We assumed no 
compression effects of the gas envelope on the solid part, a hypothe- 
sis which is justified a posteriori given the low mass fraction of gas 
obtained for each planet (see below). 

Our Bayesian model fits the planetary system as a whole, rather 
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Table 6. Parameters of the TOI-561 system, including the fifth Keplerian signal, as determined from the joint photometric and spectroscopic fit described in 


Section 5. 


Planetary parameters 


TOI-561 b TOI-561 c TOI-561 d TOI-561 e 5'^ Keplerian 

PO 0446568870099 10.77883198 000034 25102400000 770ga 473836 
To (TBJD)“ 2317.7498 + 0.0005  2238.4629+9:0008  2318.966+%:%0%3  1538.180+9-004 1664+38 
a[R. 2.685004 224392 40104535 Banu meu 
a (AU) 0.0106 + 0.0001 0.0884 + 0.0009 0.1583: 0.002 0.328 + 0.003 1.1595 
Rp/ Ry 0.0155 + 0.0004 0.0316 + 0.0004 0.0306 + 0.0008 — 0.0278*00019 - 
Rp (Ra) 1.425 + 0.037 2.91 + 0.04 2.82 + 0.07 2:55" 05 - 
b ne hee 04502 0.28013 
i (deg) 8.203 89.6901 99407021 — — S980 
Ti4 (h) 1.31 + 0.02 37502 4.543032 69804. - 
e 0 (fixed) 0.03004 34 01m. wa — uisi 
w (deg) 90 (fixed) 29153 235117 143444 348* 08 
K (ms!) 1.93 + 0.21 1.817023 3.344053 2.19 + 0.23 1.94 + 0.27 
Mp (Me) 2.00 + 0.23 Sat 13.2159 12.6 + 1.4 2043? 
Pp (Pa) 0.69 + 0.10 0.22 + 0.03 0.59 + 0.06 0.76 + 0.14 - 
Pp (g cm?) 3.84 0.5 1.2: 0.2 3.2 + 0.3 4.2+0.8 - 
Sp (Se) 4745 + 269 68.2 + 3.9 21.4 x 1.3 4.96 + 0.28 - 
T4 K) 2310 + 33 800 + 11 598 +9 415 +6 - 
gd (ms?) 9.7412 6.2+0.8 16.3 + 1.5 19.0 + 2.9 - 

Common parameters 
R,* (Ro) 0.843 + 0.005 
M,° (Mo) 0.806 + 0.036 
Px (Po) 1.31 + 0.05 
U1, TESS 0.33 + 0.08 
U2, TESS 0.23 + 0.09 
U1,CHEOPS 0.46 + 0.07 
U2,CHEOPS 0.22 + 0.09 
Harps-n (n 57) 1.401014 
rimis (ms!) 237753 
Yars n (n 5 D) 79700.41 + 0.26 
Yerrgs MS’) —1.20 + 0.42 


a TESS Barycentric Julian Date (BJD—2457000). ? Minimum mass in the hypothesis of a planetary origin. © Computed as 


Ry 
2a 


1/2 
Ta = T, (=) [f (1 — Ag)]!/^, assuming f = 1 and a null Bond albedo (Ag = 0). 4 Planetary surface gravity. © As 


determined from the stellar analysis in Section 2.1. ^ RV jitter term. € RV offset. 


than performing an independent fit for each planet, in order to ac- 
count for the correlations between the absolute planetary masses and 
radii, which depend on the stellar properties. The model fits the stel- 
lar (mass, radius, effective temperature, age, chemical abundances 
of Fe, Mg, Si), and planetary properties (RV semi-amplitudes, tran- 
sit depths, orbital periods) to derive the posterior distributions of the 
internal structure parameters. The internal structure parameters mod- 
elled for each planet are the mass fractions of the core, mantle and 
water layer, the mass of the gas envelope, the iron molar fraction in 
the core, the silicon and magnesium molar fraction in the mantle, the 
equilibrium temperature and the age of the planet (equal to the age 
of the star). For a more extensive discussion on the relation among 
input data and derived parameters we refer to Leleu et al. (2021). We 


assumed the mass fraction of the inner core, mantle, and water layer 
to be uniform on the simplex (the surface on which they add up to 
one), with the water mass fraction having an upper boundary of 0.5 
(Thiabaud et al. 2014; Marboeuf et al. 2014). For the mass of the gas 
envelope, we assumed a uniform prior in logarithmic space. Finally, 
we assumed the Si/Mg/Fe molar ratios of each planet to be equal to 
the stellar atmospheric values (even though Adibekyan et al. (2021a) 
recently showed that the stellar and planetary abundances may not 
be always correlated in a one-to-one relation). We emphasize the fact 
that, as in many Bayesian analyses, the results presented below in 
terms of planet internal structure depend to some extent on the selec- 
tion of the priors, which we chose following i.e. Dorn et al. (2017), 
Dorn et al. (2018), and Leleu et al. (2021). Analysing the same data 
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Figure 8. Phase-folded TESS (left) and CHEOPS (right) light curves of TOI- 
561 b, c, and d. Planet e shows a single transit in the TESS light curve, and it 
has no CHEOPS observations. For each planet, the coloured line indicates the 
best-fitting model, and residuals are shown in the bottom panels. Data points 
binned over 20 min (planet b) and 30 min (planets c, d and e) are shown with 
coloured dots. 


with very different priors (e.g. non uniform core/mantle/water mass 
fraction or gas fraction uniform in linear scale) would lead to different 
conclusions. 

We show the results of the internal structure modelling for the four 
planets in Figure 13. As expected from its closeness to the host star, 
planet b has basically no H/He envelope, while the other three plan- 
ets show a variable amount of gas mass. Planet c hosts a relatively 
massive gaseous envelope, with a gas mass of (5 and 95 per cent 


quantiles) Meas,c = 0.07+9:04 Me (1.303 weight percent wt%). 
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Figure 9. Phase-folded HARPS-N and HIRES RVs with residuals of TOI-561 
b, c, d and e, as resulting from the joint photometric and spectroscopic fit. 
The error bars include the jitter term added in quadrature. 


The TOI-561 system as seen by CHEOPS, HARPS-N and TESS 13 


15 
o HIRES 
T. oe — HARPS.N 
faul 
— 5 | H | 
7 | 
n 
= 0 le 
E l 
=5 4 9 
-10 
w 10 
= Oe 
E 01 et 
E 
Zo, | 
cá 8000 8700 8800 8900 9000 9100 9200 9300 9400 


BJD - 2450000 


Figure 10. Global model (grey line) with residuals of HARPS-N and HIRES 
RVs according to the 5-Keplerian photometric and spectroscopic fit. The error 
bars include the jitter term added in quadrature. 
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Figure 11. Mass-radius diagram for exoplanets with radii and masses mea- 
sured with a precision better than 30%, colour coded according to their 
incidental flux. Data are taken from the Extrasolar Planets Encyclopaedia 
catalogue (http: //exoplanet .eu/catalog/) as of 18 October 2021. The 
TOI-561 planets are labelled, and highlighted with coloured diamonds. The 
USP planets are emphasized with thick, black-contoured circles. The theoret- 
ical mass-radius curves for various chemical compositions (Zeng et al. 2019) 
are represented by solid coloured lines, while the dashed lines indicate the 
curves for an Earth-like core surrounded by a H2 envelope (2% mass frac- 
tion) at varying equilibrium temperatures. The forbidden region predicted by 
collisional stripping (Marcus et al. 2010) is marked by the shaded grey region. 
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Figure 12. Mass-radius diagram of confirmed USP planets (P < 1 d, Rp< 
2 Ra) as taken from the Extrasolar Planets Encyclopaedia catalogue in date 
18 October 2021. Iso-density lines are plotted in grey. TOI-561 b stands out 
as the lowest density USP planet known to date (pp = 3.8 + 0.5 g cm~’). 


Planet d hosts the most massive envelope (Mgas,q = 0. rus n Mo). 
which, considering the total mass of the planet, correspond to a 
smaller relative mass fraction of 0.8+10 wt%, while TOI-561 e's en- 
velope spans a range between —10.7 < log Mgas,e < —1.0, implying 
an upper limit on the gas mass of 0.11 M (< 0.9 wto). As expected 
from its low density, TOI-561 b could host a significant amount 
of water, having a water mass of My,0,b = 0.62*072 Me (81416 
wto). We stress that this result is highly dependent on the caveat 
of including only a solid water layer in the model. In fact, a mas- 
sive water layer, if present on a planet with such a high equilibrium 
temperature, would imply the presence of a massive steam atmo- 
sphere (Turbet et al. 2020). This would in turn considerably change 
the inferred water mass fraction with respect to a model that includes 
only a solid water layer. Due to the presence of the gas envelope, 
the amount of water in both planet c and d is almost unconstrained 
(My,0,¢ = 1.29*1 74 Me, i.e. 24723 wt%; Mg, oa = 3.56275 Me, 
Le. aT wto), while TOI-561 e modelling points toward a massive 


water layer, with My,0,¢ = 4.507402 Me (3615 wt%). 


7 ATMOSPHERIC EVOLUTION 


We employed the system parameters derived in this work to constrain 
the evolution of the stellar rotation period, which we use as a proxy 
for the evolution of the stellar high-energy emission affecting atmo- 
spheric escape, and the predicted initial atmospheric mass fraction of 
the detected transiting planets f*'*", that is the mass of the planetary 
atmosphere at the time of the dispersal of the protoplanetary disk. 
To this end, we used the planetary atmospheric evolution code Pasta 
described by Bonfanti et al. (2021a), which is an updated version 
of the original code presented by Kubyshkina et al. (2019b,a). The 
code models the evolution of the planetary atmospheres combining a 
model predicting planetary atmospheric escape rates based on hydro- 
dynamic simulations (this has the advantage over other commonly 
used analytical estimates to account for both XUV-driven and core- 
powered mass loss; Kubyshkina et al. 2018), a model of the stellar 
high-energy (X-ray plus extreme ultraviolet; XUV) flux evolution 
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Figure 13. Posterior distributions of the main parameters describing the internal structure of TOI-561 b (top left), c (top right), d (bottom left), and e (bottom 
right). Each corner plot shows the mass fraction of the inner core and of the water layer, the molar fractions of silicon and magnesium in the mantle, the iron 
molar fraction in the inner core, and the mass of gas in logarithmic scale. On top of each column are printed the mean and the 5 per cent and 95 per cent quantiles 
values. For each planet, and we show an illustration of the radius fractions of the inner core+mantle (dark gray), water layer (dark blue), and gas envelope (light 
blue), corresponding to the medians of the posterior distributions. The coloured rectangles indicate the uncertainty on the corresponding layer thickness, while 
the black dashed outer rings represent the uncertainty on the total radius. Equilibrium temperature and planetary surface gravity are reported for each planet. 


(Bonfanti et al. 2021a), a model relating planetary parameters and toplanetary disk, and that the planets hosted at some point in the past 
atmospheric mass (Johnstone et al. 2015b), and stellar evolutionary or still host a hydrogen-dominated atmosphere. 

tracks (Choi et al. 2016). The main assumptions of the framework 

are that planet migration did not occur after the dispersal of the pro- For each planet, the evolution calculations begin at an age of 


5 Myr, which is the age assumed in the code for the dispersal of 
the protoplanetary disk. At each time step, the framework derives 


MNRAS 000, 1-22 (2022) 


The TOI-561 system as seen by CHEOPS, HARPS-N and TESS 15 


the mass-loss rate from the atmospheric escape model employing 
the stellar flux and the system parameters, and uses it to update 
the atmospheric mass fraction. This procedure is then repeated until 
the age of the system is reached or the planetary atmosphere has 
completely escaped. The free parameters of the algorithm are the 
initial atmospheric mass fraction at the time of the dispersal of the 
protoplanetary disk, and the indexes of the power law controlling 
the stellar rotation period (see Bonfanti et al. 2021a, for a detailed 
description of the mathematical formulation of the power law), that 
we use as proxy for the stellar XUV emission. 

The free parameters are constrained by implementing the atmo- 
spheric evolution algorithm in a Bayesian framework employing the 
MCMC tool presented by Cubillos et al. (2017). The framework uses 
the system parameters with their uncertainties as input priors. It then 
computes millions of forward planetary evolutionary tracks, varying 
the input parameters according to the shape of the prior distributions, 
and varying the free parameters within pre-defined ranges, fitting the 
current planetary atmospheric mass fractions obtained as described 
in Section 6. The fit is done at the same time for all planets, thus 
simultaneously constraining the rotational period, and the results are 
posterior distributions of the free parameters. In particular, we opted 
for fitting for the planetary atmospheric mass fractions instead of the 
planetary radii. This enables the code to be more accurate by avoid- 
ing the continuous conversion of the atmospheric mass fraction into 
planetary radius, given the other system parameters (see also Delrez 
et al. 2021). 

Figure 14 shows the results of the planetary atmospheric evolu- 
tion simulations. As a proxy for the evolution of the stellar rotation 
period, in Figure 14, we show the posterior distribution of the stellar 
rotation period at an age of 150 Myr, further comparing it to the 
distribution of stellar rotation periods observed in stars member of 
young clusters of comparable age and with masses that deviate from 
M. less than 0.1 Mo (from Johnstone et al. 2015a). The inferred 
posterior distribution for the rotation period is consistent with mem- 
bership of the slowly-rotating period-colours sequence in clusters of 
this age. However, this comparison should be taken with some cau- 
tion, since there are no comprehensive studies on the rotation-colour 
distributions of 150 Myr-old clusters with the same metallicity as 
TOI-561. The initial atmospheric mass fractions of planets b and c 
are rather broad and peak at about one planetary mass. This is be- 
cause both planets are close enough to the host star and have a small 
enough mass to have been subject to significant atmospheric escape. 
Therefore, to enable the presence of a thin hydrogen atmosphere, as 
predicted by the internal structure model, both planets had to host a 
significant hydrogen envelope after the formation and atmospheric 
accretion processes. Instead, planets d and e are far from the host 
star and massive enough not to have been subject to significant atmo- 
spheric escape, which is why we obtain an initial atmospheric mass 
fraction that resembles the current one. We also find that the posterior 
distributions of all input parameters match well the inserted priors 
(not shown here). As a whole, the results indicate that the currently 
Observed system parameters are compatible with a scenario in which 
migration happened (if at all) exclusively inside the protoplanetary 
disk. Otherwise the code would have led to mismatches between 
the prior and posterior of the input parameters (particularly for what 
concerns the planetary masses and/or the stellar mass and age), in 
addition to showing incoherent results in the posterior distribution of 
the output parameters. This is for example the case of the TOI-1064 
system, which is composed by two transiting planets with comparable 
masses and irradiation levels, but significantly different radii (Wilson 
et al. 2022). In our framework in which planets do not migrate after 
the dispersal of the protoplanetary nebula, reproducing the physical 


parameters of the planets composing the TOI-1064 system requires 
different evolutions of the stellar rotation rate, which is not possible, 
thus calling for a post-nebula migration. 


8 DISCUSSION AND CONCLUSIONS 


In this study, we confirm the presence of four transiting planets 
around TOI-561, with orbital periods of approximately 0.44, 10.8, 
25.7, and 77 days (Table 6). Our analysis disproves the presence of the 
previously suggested planet TOI-561 f (P ~ 16.3 day; W21). TOI- 
561 is one of the few 4-planet systems having precise radius and mass 
measurements for all the planets. Thanks to our global photometric 
and RV analysis, we refined all masses and radii with respect to the 
L21 values, and we precisely determined the planetary bulk densities, 
with uncertainties of 14.4%, 13.6%, 10.2%, and 18.4% for planets b, 
c, d, ande, respectively. The higher uncertainty on planet e reflects the 
lower precision in the radius determination (5% uncertainty), which 
is based on the analysis of a single TESS transit, and highlights the 
importance of the high-precision CHEOPS photometry. In fact, with 
asingle CHEOPS transit we managed to decrease the uncertainty on 
the radius of planet d from 5.1% (L21, based on one TESS transit) 
to 2.5%. Including also the improvement on the mass, this implied a 
decrease on the density uncertainty from 18.9% to 10.2%. We expect 
asimilar improvement for planet e with future CHEOPS observations 
scheduled for 2022. The improvement in the radius of TOI-561 e is 
particularly important, since the planet is an interesting target for the 
study of the internal structure of cold sub-Neptunes. Its long period 
(Pa = 71.033025 d) implies an insolation flux of Se = 4.96:-0.28 Se 
and a relatively cool zero Bond albedo equilibrium temperature of 
Teg,e = 415 + 6 K. As shown in Figure 15, TOI-561 e is one of the 
few cool, long-period planets orbiting a star bright enough for precise 
RV characterization, and it is therefore an optimal test-case to refine 
tools and models that will be useful to characterize targets of future 
long-staring missions like PLATO. 

TOI-561 hosts one of the most intriguing USP planets discov- 
ered to date. As initially suggested by L21, our analysis confirms 
that TOI-561 b is the lowest density (pp = 3.8 + 0.5 g cm?) USP 
super-Earth that we know of (see Figure 12), and it paves the way for 
in-depth studies of interior composition, and formation and evolu- 
tion processes of USP planets. Even though now the mass values are 
consistent within lo, contrary to what proposed by W21 (see Sec- 
tion 2.2) TOI-561 b is not consistent with a pure rocky composition, 
and to explain the planetary density our internal structure modelling 
(Section 6) predicts basically no H/He envelope, and a massive water 
layer. In this regard, an important point to consider is that, with an 
insolation flux of Sy ~ 4745 Sẹ, the planet receives more irradiation 
from the star than the theoretical runaway greenhouse limit (Kasting 
et al. 1993; Goldblatt & Watson 2012; Kopparapu et al. 2013). In this 
case, a large water content would imply the presence of an extended 
steam atmosphere, which in turn would increase the measured ra- 
dius with respect to a purely condensed water world, leading in our 
model to an overestimation of the bulk water content (Turbet et al. 
2020). The presence of a water steam envelope could eventually be 
tested with the James Webb Space Telescope (JWST). In fact, with an 
Emission Spectroscopy Metric (ESM, Kempton et al. 2018) value of 
8.2, TOI-561 b is a promising target for secondary eclipse and phase 
curve observations. More complex models, including a lighter core 
compositions (i.e. a Ca/Al enriched core), the modelling of water 
steam envelopes, or wet-melt solid interiors related to deep water 
reservoirs (Dorn & Lichtenberg 2021), could be an interesting step 
forward in the understanding of the planet structure and composi- 
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Figure 14. From left to right: Posterior distributions (dark blue lines) of the stellar rotation period at an age of 150 Myr and of the initial atmospheric mass 
fractions of TOI-561 b, c, d and e. In each panel, the purple region represents the 68 per cent highest probability density intervals. In the left panel, the black 
thin line shows the rotation period distribution of stars member of open clusters with ages around 150 Myr. Data are taken from Johnstone et al. (2015a), who 
report the rotation period of ~ 2000 stars belonging to the Pleiades, M50, M35, and NGC 2516, whose ages are between 125 and 150 Myr. To generate the 
black histogram we selected a sub-sample of 578 stars, which have masses that deviate from M, less than 0.1 Mo. In the other panels, the horizontal orange 
lines mark the uniform prior used in the fit, scaled to the highest peak of each posterior distribution for better visualization. The light blue lines indicate the 
current atmospheric mass fraction of each planet determined as described in Section 6. 
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Figure 15. V magnitude versus planetary periods for confirmed transiting 
exoplanets as reported in the Extrasolar Planets Encyclopaedia catalogue 
(http://exoplanet.eu/catalog/)in date 18 October 2021. The dashed 
vertical line marks V — 12 mag. TOI-561 e is one of the few long-period 
planets orbiting a star bright enough for precise RV characterization. 


tion. The low density of TOI-561 b could also be related to the fact 
that the host star is a metal-poor, thick-disk star. Adibekyan et al. 
(2021a) showed that the composition of the rocky planets reflects the 
chemical abundances of the host star (even though not in a one-to- 
one relation), so implying a lighter composition for TOI-561 b with 
respect to other USP planets that orbit more metal-rich stars!?. Ac- 
cording to Adibekyan et al. (20212), the low density of TOI-561 b is 
consistent with the general p/PẸarth-like — Fed trend and dispersion 
inferred from the sample of rocky planets analysed by the authors 
(see Figs. 2, 3 therein), where p/PẸarth-like 1s the planetary density 


normalised to that expected for an Earth-like composition, and i 


13 All the USP planets shown in Figure 11 have [Fe/H]> —0.14. 
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is the iron-to-silicate mass fraction of the protoplanetary disk as in- 
ferred from the stellar properties. An additional interesting remark 
concerns the Galactic kinematics of the host star. According to our 
analysis, performed as described in Mustill et al. (2021), TOI-561 is 
located in a low-density region of the 6-dimensional Galactic phase 
space (see Winter et al. 2020, Mustill et al. 2021, and Kruijssen 
et al. 2021 for definition and discussion), which is not surprising 
given that TOI-561 is a thick disk star (Mustill et al. 2021). Kruijssen 
et al. (2020) showed that stars in low-density regions seem to host 
no super-Earths, but only sub-Neptunes, i.e. planets having a signif- 
icant H/He envelope and therefore located above the radius gap. In 
this context, TOI-561 b is an interesting object that runs counter to 
this finding. We point out that this result should be taken with some 
caution, since the Kruijssen et al. (2020) sample does not include 
planets with periods shorter than one day, and it excludes stars with 
ages » 4.5 Gyr'^. 

All the four planets seem to host a large water layer (Section 6), 
although with high uncertainties, especially for planet c and d, due 
to the degeneracy related to the possible presence of a gas envelope. 
Also in this case, the presence of a considerable amount of water 
could be linked with the stellar properties. In fact, Santos et al. 
(2017) showed that metal-poor, thick disk stars are expected to form 
planetary building blocks with a higher water mass fraction (~ 76%) 
compared to metal-rich, thin disk stars (^ 58%). Therefore, we would 
expect these stars to produce water-rich planets, a result that is in 
agreement with our findings on the TOI-561 system. 

Except for TOI-561 b, all the other planets are suggested to host a 
non-negligible H/He envelope. In particular, the gas content of planet 
€ (~ 1.3 wt%, the highest mass fraction among the four planets) im- 
plies a much lower density with respect to the density of planet d, 
even though the two planets have a similar size. This is reflected 
in the different positions of the planets in the mass-radius diagram 
(Figure 11). The two planets show hints of a different evolution for 
what concerns their gas content. In fact, our atmospheric evolution 
model (Section 7) suggests that planet c underwent a strong enve- 
lope loss after the atmospheric accretion and the dispersal of the 
protoplanetary nebula, while planet d (as well as planet e) did not 


1^ We note however that the stellar ages used in Kruijssen et al. (2020) are 
quite inhomogeneous, coming directly from the NASA Exoplanet Archive, 
and can therefore show a large scatter with respect to a homogeneous deter- 
mination (Adibekyan et al. 2021b). 
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experience strong atmospheric escape, with a current gas content that 
is comparable to the original one. The surprising difference in gas 
mass fraction between planets c, d and e, not only at present time but 
also at the end of their formation phase, takes probably its origin in 
the conditions that prevailed during the protoplanetary disk phase. 
Planet c is indeed likely sub-critical because of its low mass, where 
sub-critical planets are those with masses below the critical value 
required to initiate runaway gas accretion (see Helled et al. 2014 for 
a recent review on the core accretion model), whereas planets d and 
e never accreted large amounts of gas as demonstrated in Section 7, 
and so they also remained always below the critical mass. The in- 
terpretation of the different gas mass fractions could therefore result 
from the structure of sub-critical planets. In this case, the gas mass 
fraction depends on the core mass, the thermodynamical properties 
in the disk, and more importantly the accretion rate of solids (lower 
accretion rate translating in larger gas mass fraction). Interpreting 
the internal structure of the four planets of the system in a global 
planetary system formation model could therefore constrain these 
parameters. 

With its derived properties, TOI-561 c has a Transmission Spec- 
troscopy Metric (TSM, Kempton et al. 2018) of 110.4, and is there- 
fore a suitable target for atmospheric characterization with J WST.^ 
Instead, planets d and e have lower TSM values of 30.7 and 16.2, 
respectively. As the TSM is proportional to the equilibrium temper- 
ature, it is not surprising to obtain lower values for the two planets, 
given their longer periods. 

In addition to the characterization of the four planets, we also 
identified a significant long-term signal (P ~ 473 d) in the RVs. 
On the basis of our current dataset, we are not able to distinguish 
between a stellar (magnetically-induced) or planetary origin. Long- 
term monitoring using both spectroscopic ground-based facilities 
and future long-staring missions like the PLATO spacecraft will 
allow us to shed light on the nature of this additional signal, and 
to potentially find new outer companions. It is worth noting that, 
if the above-mentioned signal proves in future to be of planetary 
origin, there is a non-zero chance that, under the assumption of co- 
planarity, such a planet would transit. In fact, assuming the same 


inclination of planet e and using the semi-major axis a/ Ry — pe Lees 


derived from our global fit, we infer an impact parameter of UNE 


Moreover, the planet would orbit in TOI-561's empirical habitable 
zone (175 x P sz 652 d), as originally defined by Kasting et al. 
(1993) using a 1D climate model, and later updated in Kopparapu 
et al. (2013); Ramirez & Kaltenegger (2016) for main-sequence stars 
with 2600 < Tog < 10000 K. 

This work bears witness to the fruitful results that can be obtained 
by the timely combination of data coming from different instruments. 
It adds to the works (Bonfanti et al. 2021b; Leleu et al. 2021; Delrez 
et al. 2021) that prove the potential of CHEOPS in precisely char- 
acterizing TESS-discovered exoplanets, as well as demonstrating the 
key role of high-precision spectrographs such as HARPS-N when 
working in synergy with space-based facilities. 
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APPENDIX A: CHEOPS LIGHT CURVES AND 
TELEGRAPHIC PIXEL TREATMENT 


As described in Section 3.2, the three CHEOPS visits of TOI-561 
were reduced via the standard DRP processing. The light curves 
presented in this study, obtained using the RINF aperture size (RINF 
= 0.9 x DEFAULT, where DEFAULT = 25 px; see also Section 3.2), 
are shown in Figure A1. While for the two initial visits the automatic 
DRP processes was performed, the appearance of some telegraphic 
pixels during the third visit required a more in-depth analysis. 

In addition to the large number of known hot pixels present in 
the CHEOPS CCD (some of them visible in Figure 2), some nor- 
mal pixels can change their behaviour during the duration of a visit, 
for example becoming ‘hot’ after a SAA crossing. These pixels, 
called ‘telegraphic’ for their abnormal behaviour, can affect the pho- 
tometry if located within the photometric aperture (see for exam- 
ple Leleu et al. 2021). During the third CHEOPS visit, we identi- 
fied an unusual flux bump before the ingress of TOI-561 d transit, 
at BJD ~ 2459318.75 (top panel, Figure A2). After analyzing the 
statistics of each pixel light curve within the photometric aperture, 
we detected a telegraphic pixel with a large flux variation (second 
panel, Figure A2) located within the CHEOPS PSF. The exact po- 
sition of this pixel on the CHEOPS CCD is shown in Figure 2. We 
masked the pixel flux and repeated the photometric extraction of the 
visit using the RINF aperture, so removing the flux jump in the light 
curve (bottom panel, Figure A2). During this analysis, we detected 
two additional telegraphic pixels within the photometric aperture, in- 
ducing smaller, but still significant variations in the light curve flux 
(third panel, Figure A2). We corrected for the effect of these pixels 
as described above. 

While investigating the nature of the flux bump happening dur- 
ing the third visit, we also extracted the light curve using a PSF- 
photometry approach exploiting the PIPE (PSF Imagette Photomet- 
ric Extraction) software!Ó. PIPE is a photometric extraction package 
specifically developed to extract CHEOPS light curves by applying 
PSF photometry on the 60-pixel imagettes, complementing the offi- 
cal DRP extraction. The use of PSF photometry makes usually easier 
to filter out the impact of hot pixels and cosmic rays, by either giving 
them a lower weight or masking them entirely in the fitting pro- 
cess. However, in this case the telegraphic pixel was located inside 
the CHEOPS PSF, requiring a careful manual masking. As for the 
DRP light curve, the flux bump in the PIPE photometry is reduced 
after masking the telegraphic pixel (bottom panel, Figure A3). The 


16 https://pipe-cheops.readthedocs.io/en/latest/index.html 
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Figure A1. CHEOPS RINF light curves of TOI-561 as extracted from the 
DRP, with 40 -clipping for outliers removal. Visits 1, 2 and 3 are shown from 
top to bottom. 


PIPE-extracted light curve resulted in a slightly lower mean absolute 
deviation (MAD) with respect to the DRP photometry (top panel, 
Figure A3), mainly due to the lower number of outliers present in the 
PSF photometry. For a more detailed comparison between PIPE and 
DRP photometries, see Morris et al. (2021b). We performed the same 
global analysis described in Section 5 using the PIPE light curve in- 
stead of the DRP one, obtaining consistent results and comparable 
uncertainties on the transit parameters of both planets b and d. We 
therefore decided to use the light curve obtained with the official 
DRP extraction in our final analysis. 


APPENDIX B: FLOATING CHUNK OFFSET METHOD ON 
TOI-561 B 


In order to investigate the literature discrepancy on the mass of 
TOI-561 b (Section 2.2), we adopted a specific observing strategy 
with HARPS-N targeting the USP planet (Section 3.3), obtaining 
multiple observations during the same night for 22 nights. Multiple 
nightly observations can be used to precisely infer the mass of USP 
planets using the Floating Chunk Offset method (FCO; Hatzes 2014), 
which consists in applying a nightly offset to remove all the other 
signals present in the system, both of planetary and stellar origin (i.e. 
Howard et al. 2013; Pepe et al. 2013; Malavolta et al. 2018; Frustagli 
et al. 2020). The FCO method is only applicable when the separation 
between the USP period and the period of all the other signals is large 
enough, and the RV semi-amplitude has a similar or larger value with 
respect to the other signals. As demonstrated in L21, these conditions 
apply to TOI-561 b, for which the authors derived an FCO semi- 
amplitude of Ky pco = 1.8020.38 m s7! (My Fco = 1.83+0.39 Ma) 
exploiting multiple observations collected over ten nights. 

Here, we applied the FCO method to TOI-561 b on a total of 22 
HARPS-N nights, adding 12 novel nights to the 10 nights already 
presented in L21. Out of the total set, four nights have six multiple 
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Figure A2. Top panel: TOI-561 RINF original light curve of the third visit 
(light blue dots) after the removal of 4c outliers, with over-plotted the 15- 
minute binned light curve (dark blue dots). The start of the flux jump due to 
the telegraphic pixel is marked with the dashed vertical line. Second panel: 
light curve of the telegraphic pixel located within the CHEOPS PSF. Third 
panel: light curve of the two additional telegraphic pixels located within the 
RINF aperture. Bottom panel: corrected light curve after masking the three 
telegraphic pixels. 


Observations extending over more than 40 per cent of the orbital 
period of the planet, and span opposite orbital phases to provide 
an optimal phase coverage. We performed a PyDE + emcee fit with 
PyORBIT, assuming a fixed zero eccentricity and Gaussian priors 
on period and Tọ coming from the global fit, and we added a jit- 
ter term to account for possible additional white noise. We derived 
a semi-amplitude of Ky = 1.81 + 0.31 m st, corresponding to a 
mass of My = 1.86 + 0.33 Me, with a jitter of 0.96902 m sl. 
Figure B1 shows the resulting phase-folded RVs. The derived mass 
and semi-amplitude are nicely in agreement with the L21 values, 
and they support the values inferred from our joint photometric and 
RV modelling (Section 5), being consistent within lo. Given the 
higher number of RVs included in the joint fit, which led to smaller 
uncertainties on the derived parameters, we decided to adopt as final 
values for TOI-561 b the ones obtained from the global modelling, 
ie. Kp = 1.93 + 0.21 ms^!, My = 1.99 + 0.22 Mo. 


1.0025 


« 1.0000 


0.9975 


Relative flux 


0.9950 


0.9925 


95 


X 
y 
2 
c 
N 
[S 


PIPE corrected 
ET 


£o 


eire alatus! 


1.0000 


elative flus 


R 
o 
io 
e 
x 
e 


0.5 1.0 15 2.0 25 
BJD-2459317 


Figure A3. Top panel: comparison between DRP and PIPE-extracted light 
curve of TOI-561 third visit, before the telegraphic pixel correction and with 
4c outliers removal. The DRP has an MAD of 371 ppm over the whole visit, 
while PIPE of 325 ppm. Bottom panel: PIPE light curve after the telegraphic 
pixel correction. The light curve gets slightly noisier (MAD = 331 ppm) 
because one less pixel is considered in the reduction, but more reliable thanks 
to the exclusion of the telegraphic pixel flux. In both panels, the beginning of 
the flux jump is highlighted with a vertical dashed line. 
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Figure B1. Phase-folded RVs of the 22 HARPS-N nights used to model TOI- 
561 b with the FCO method. The error bars include the jitter term added in 
quadrature. 
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